A fracture model with variable range of interaction 
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We introduce a fiber bundle model where the interaction among fibers is modeled by an adjustable 
stress-transfer function which can interpolate between the two limiting cases of load redistribution, 
the global and the local load sharing schemes. By varying the range of interaction several features 
of the model are numerically studied and a crossover from mean field to short range behavior is 
obtained. The properties of the two regimes and the emergence of the crossover in between are 
explored by numerically studying the dependence of the ultimate strength of the material on the 
system size, the distribution of avalanches of breakings, and of the cluster sizes of broken fibers. 
Finally, we analyze the moments of the cluster size distributions to accurately determine the value 
at which the crossover is observed. 
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I. INTRODUCTION 

Fracture processes have attracted the attention of the 
scientific community since many years. Processes involv- 
ing heterogeneous systems, for which a definite and com- 
plete physical description has not been found despite the 
many partial successes of the last decades, are of special 
theoretical and practical interest. In particular, 

the latest developments of statistical mechanics have led 
to a deeper understanding of breakdown phenomena in 
heterogeneous systems, but some fundamental questions 
remain unsolved. The difficulties arise because in model- 
ing fracture of heterogeneous materials, one has to deal 
with systems formed by many interacting constituents, 
each one having different statistical properties related to 
some breaking characteristics of the material, distributed 
randomly in space and/or time EE]- So, the complete 
analytical solution is in almost all cases prohibitive and 
one has to solve the problem by means of numerical sim- 
ulations or to study simplified models which can be an- 
alytically tractable (at least in some limits) in order to 
gain physical insights that guide our understanding to 
more complex models. 

The major challenge in dealing with fracture problems 
is to combine the statistical evolution of damage across 
the entire macroscopic system and the associated stress 
redistributions to accurately predict the point of final 
rupture of the material. In doing this linkage, one has 
to take care in order not to make major simplifications 
particularly in the redistribution rule, where a great deal 
of the physics of the problem is hidden. A very impor- 
tant class of approaches to the fracture problem are the 
well-known Fiber Bundle Models (FBM), which were in- 
troduced long time ago by Daniels [Q and Coleman ||] 
and have been the subject of intense research during the 
last several years FBM's are constructed so that 



a set of fibers is arranged in parallel each one having a 
statistically distributed strength. The specimen is loaded 
parallel to the fiber direction and the fibers break if the 
load acting on them exceeds their threshold value. Once 
the fibers begin to fail one can choose among several load 
transfer rules. The simplest case is to assume global load 
sharing (GLS) which means that after each fiber failure, 
its load is equally redistributed among all the intact fibers 
remaining in the set. This model, known as global fiber 
bundle model, is a mean field approximation where long 
range interactions among the elements of the system are 
assumed and can be solved analytically [ pdp8| . At the 
other extreme, one finds the local load sharing{LLS) fiber 
bundle model where the load borne by failing elements 
is transferred to their nearest neighbors. This case rep- 
resents short range interactions among the fibers. Other 
schemes have been proposed with a relative success in de- 
scribing rupture processes at large scale like earthquakes 
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Despite their simplicity, FBM are very important be- 
cause they capture most of the main aspects of material 
damage and breakdown. They have provided a deeper 
understanding of fracture processes and have served as 
a starting point for more complex models of fiber re- 
inforced composites and other micro-mechanical models 
p2|-^|. However, stress redistribution in actual hetero- 
geneous materials should fall somewhere in between LLS 
and GLS since there is an important fraction of stress re- 
distributed to other intact elements not localized in the 
neighborhood of the failed ones, nevertheless maintain- 
ing stress concentrations around the broken fibers. With 
this aim, several studies have been carried out during the 
last two decades and Monte Carlo simulations have been 
used to numerically study the distribution of composite 
strengths in 2D and 3D for different fiber arrangements 
Jf6| , p5| -[2S|] . Nevertheless, in order to obtain reliable con- 
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elusions the number of fibers forming the system has to be 
very large which makes the numerical problem, in many 
cases, too time consuming as to perform the study in a 
reasonable amount of time. 

In this paper, we introduce a fiber bundle model where 
the interaction among fibers is modeled by an adjustable 
stress-transfer function, which interpolates between the 
two limiting cases of load redistribution, the global and 
the local load sharing schemes. By varying the effective 
range of interaction one observes a crossover from mean 
field to short range behavior. To explore the properties 
of the two regimes and the emergence of the crossover in 
between, a comprehensive numerical study of the model 
is performed. We study the dependence of the ultimate 
strength of the material on the system size and found 
that the system has only one nonzero critical load in the 
thermodynamic limit. When no critical point exits, the 
ultimate strength of the material goes to zero exactly as 
in local load sharing models as ~ ln ^ N ^ , with increas- 
ing system size N. We also study the distribution of 
avalanches of fiber breaks, and that of the cluster sizes 
of broken fibers for the two distinct regimes and perform 
a moment analysis to accurately determine the crossover 
value. 

The rest of the paper is organized as follows. The next 
section is devoted to present the model and the way in 
which numerical simulations are carried out. The results 
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obtained are presented and analyzed in Section III. The 
final section is devoted to discussions and to state our 
conclusions. 



II. THE MODEL 

The fracture of heterogeneous systems is character- 
ized by the highly localized concentration of stresses at 
the crack tips that makes possible the nucleation of new 
cracks at these regions such that the actual crack grows 
leading to the final collapse of the system. In elastic ma- 
terials, the stress redistribution follows a power law, 
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where <r a dd is the stress increase on a material element 
at a distance r from the crack tip. The above general 
relation covers the cases of global and local load sharing, 
widely used in fiber bundle models of fracture, as the 
limiting cases 7 — ► 0, and 7 — > 00, respectively. 

Motivated by the above result of fracture mechanics 
we introduce a fiber bundle model where the load shar- 
ing rule takes the form of Eq. (|l|). Suppose a set of N 
parallel fibers each one having statistically distributed 
strength taken from a probability distribution function 
P and identified by an integer i, 1 < i < N. In materials 
science, the Weibull distribution has been proved to be 
a good empirical statistical distribution for representing 
fiber strength, 



where p is the so-called Weibull index, which controls the 
degree of threshold disorder in the system (the bigger 
the Weibull index, the narrower the range of threshold 
values), and ao is a reference load which acts as unity. 
Thus, to each fiber i a random threshold value cr ? - is 
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assigned. The system is driven by increasing quasistati- 
cally the load on it, which is performed by locating the 
fiber which minimizes <Ji — <Ji th and adding this amount of 
load to all the intact fibers in the system. This provokes 
the failure of at least one fiber which transfers its load 
to the surviving elements of the set. This may provoke 
other fractures in the system which in turn induce ter- 
tiary ruptures and so on until the system fails or reaches 
an equilibrium state where the load on the intact fibers 
is lower than their individual strengths. In this later 
case, the slow external driving is applied again and the 
process is repeated up to the macroscopic failure of the 
material. The number of broken fibers between two suc- 
cessive external drivings is the size of an avalanche s, and 
the number of parallel updatings of the lattice during an 
avalanche is called its lifetime T. 

We now focus on the load transfer process following 
fiber failures. We suppose that, in general, all intact 
fibers have a nonzero probability of being affected by the 
ongoing failure event, and that the additional load re- 
ceived by an intact fiber i depends on its distance 
from fiber j which has just been broken. Furthermore, 
elastic interaction is assumed between fibers such that 
the load received by a fiber follows the power law form of 
Eq. (|l|). Hence, in our discrete model the stress-transfer 
function ^(7*^,7) takes the form 
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where 7 is our adjustable parameter, Z is given by the 
normalization condition Z = (Yl r~^)~ l (the sum runs 
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over the set / of all intact elements) and r*y is the 
distance of fiber i to the rupture point (xj,yj), i.e., 
r-ij = -J (xi — Xj) 2 + (yi — yf) 2 in 2D. Periodic bound- 
ary conditions are assumed so that the largest r value is 

R-max = ^o" 1 ^ ' wnere L 1S the linear size of the system. 
We note here that the assumption of periodic boundary 
conditions is for simplicity. In principle, an Ewald sum- 
mation procedure would be more accurate. The model 
construction is illustrated in Fig. [j]. It is easy to see that 
in the limits 7 — > and 7 — > 00 we recover the two ex- 
treme cases of load redistribution in fiber bundle models: 
the global load sharing and the local load sharing, re- 
spectively. We should note here that, strictly speaking, 
for all 7 different from the two limits above, the range 
of interaction covers the whole lattice. However, when 
changing this exponent, one moves from a very localized 
effective range of interaction to a truly global one as 7 
approaches zero. So, we will refer henceforth to a change 
in the effective range of interaction. 
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In summary, during an avalanche of failure events, an 
intact fiber i receives at each time step r the load borne 
by failing elements j. Consequently, its load increases by 
an amount, 

o-i(t + r) =ai(t + T-l)+ °'j{ t + T - 1 ) F i r ihl)^ 

j£B(r) 

(3) 

where the sum runs over the set B(t) of elements that 
have failed in a time step r. Thus, <Ji(to + T) = 

T 

Y] Ui (to + t) is the total load element i receives during 

r=l 

an avalanche initiated at to and which ended at t a + T. 
In this way, when an avalanche ends, the external field 
is applied again and another avalanche is initiated. The 
process is repeated until no intact elements remain in 
the system and the ultimate strength of the material a c , 
is defined as the maximum load the system can support 
before its complete breakdown. 

Unfortunately, the complete analytical approach to the 
general model introduced here is inaccessible. There are 
a few cases where this task can be achieved such as the 
global load sharing model where the load acting on sur- 
viving elements for a given external force F is known 
jl4],[l8| and some ID models [f30|-[32[ (which are irrelevant 
for practical purposes) . The main difficulty is that in or- 
der to analytically solve the problem, one needs to know 
the transition probabilities for all the possible paths lead- 
ing the system from the state in which all the elements 
are intact to the state in which they have failed. This cal- 
culation eventually becomes impossible for large system 
sizes. So, a first step is to learn from Monte Carlo simula- 
tions which, furthermore, allows us to better understand 
the physical mechanisms of fracture and to study mod- 
els difficult to handle analytically as well as to guide our 
search for analytical calculations. 



III. MONTE CARLO SIMULATION OF THE 
FAILURE PROCESS 

We have carried out large scale numerical simulations 
of the model described above in two dimensions. The 
fibers are identified with the sites of a square lattice of 
linear size L with periodic boundary conditions. The 
failure process is then simulated by varying the effective 
range of interaction between fibers by controlling 7, and 
recording the avalanche size distribution, the cluster size 
distribution and the ultimate strength of the material 
for several system sizes. Each numerical simulation was 
performed over at least 50 different realizations of the 
disorder distribution. 

Figure |2| shows the ultimate strength of the material 
for different values of the parameter 7 and several sys- 
tem sizes from L = 33 to L — 257. Clearly, two dis- 
tinct regions can be distinguished. For small 7, a c is 



independent, within statistical errors, of both the effec- 
tive range of interaction and the system size. At a given 
point 7 = 7 C a crossover is observed, where 7 C falls in the 
vicinity of 7 = 2. The region 7 > j c might eventually be 
further divided into two parts, the first region character- 
ized by the dependence of the ultimate strength of the 
bundle on both the system size and the effective range of 
interaction; and a second region where a c only depends 
on the system size. This would mean that there might be 
two transition points in the model, for which the system 
displays qualitatively and quantitatively different behav- 
iors. For 7 < 7 C the ultimate strength of the bundle 
behaves as in the limiting case of global load sharing, 
whereas for 7 > 7 C the local load sharing behavior seems 
to prevail. Nevertheless, the most important feature is 
that when decreasing the effective range of interaction in 
the thermodynamic limit, for 7 > j c , the critical load 
is zero. This observation is further supported by Fig. |^, 
where we have plotted the evolution of a c as a function 
of for different values of the exponent 7. Here, the 
two limiting cases are again clearly differentiated. For 
large 7 all curves decreases when N — > 00 as 
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This qualifies for a genuine short range behavior as found 
in LLS models where the same relation was obtained for 
the asymptotic strength of the bundle Jl^,^0| . It is worth 
noting that in the model we are analyzing, the limit- 
ing case of local load sharing corresponds to models in 
which short range interactions are considered to affect 
the nearest and the next-nearest neighbors. In the tran- 
sition region, the maximum load the system can support 
also decreases as we approach the thermodynamic limit, 
but in this case much slower than for 7 S> 7 C - It has been 
pointed out that for some modalities of stress transfer, 
which can be considered as intermediate between GLS 
and LLS, a c decreases for large system sizes following 
the relation a c ~ i n (iniv) as m ^ e case °f hierarchical 

. In our case, we have fitted 



ln(lniV) 

load transfer models |33j 
our results with this relation but we have not obtained 
a single collapsed curve because the slopes continuously 
vary until the LLS limit is reached. Finally, the region 
where the ultimate stress does not depend on the system 
size shows the behavior expected for the standard GLS 
model, where the critical load can be exactly computed 
as a c — (pe)~ 1 / p for the Weibull distribution. The nu- 
merical values obtained for p = 2 are in good agreement 
with this later expression. 

The fracture process can also be investigated by look- 
ing at the recursory activity before the complete break- 
down. The statistical properties of rupture sequences are 
characterized by the avalanche size distribution which 
from the experimental point of view could be related to 
the acoustic emissions generated during the fracture of 
materials p^-p7| . Figure || shows the avalanche size dis- 
tribution for different values of 7. Again, we observe that 
for decreasing effective range of interaction (increasing 7) 
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there is a crossover in the distribution of avalanche sizes. 
The upper curves can be very well fitted by a power law 
P{m) ~ m~ T , with r ~ §, the value obtained for long 
range interactions ||l7|-^o||. 

As soon as the localized nature of the interaction be- 
comes dominant 7 > j c , the power law dependence of 
the avalanche size distribution with the exponent t « § 
does not apply anymore. The lack of a characteristic 
size is a fingerprint of a highly fluctuating activity that 
could be related to the very nature of the long range 
interactions. The avalanche size distribution is a mea- 
sure of causally connected broken sites and the spatial 
correlations in this limit are ruled out. All the intact ele- 
ments have a nonzero chance to fail independently of the 
(spatial) rupture history, and any given element could 
be near to its rupture point regardless of its position in 
the lattice. This is not indeed the case when 7 is large 
enough and the short range interaction prevails. Now, 
the spatial correlations are important and concentration 
of stress takes place in the fibers located at the perimeter 
of an already formed cluster. Fibers far away the clus- 
ters of broken elements have significantly lower stresses 
and thus the size of the largest avalanche is reduced as 
well as the number of failed fibers belonging to the same 
avalanche, leading to a lower precursory activity. 

A further characterization of what is going on in the 
fracture process can be carried out by focusing on the 
properties of clusters of broken fibers. The clusters 
formed during the evolution of the fracture process are 
sets of spatially connected broken sites on the square lat- 
tice ||[l2]]. We consider the clusters just before the global 
failure and they are defined taking into account solely 
nearest neighbors connections. It is important to note 
that the case of global load sharing does not assume any 
spatial structure of fibers since it corresponds to the mean 
field approach. However, in our case it is obtained as a 
limiting case of a local load sharing model on a square 
lattice, which justifies the cluster analysis also for GLS. 
Fig. H illustrates how the cluster structure just before 
complete breakdown changes for various values of 7. We 
have also recorded the cluster size distribution as a func- 
tion of the effective range of interaction. Figure |^ shows 
the size distribution n(s c ), of the two-dimensional clus- 
ters for several values of the exponent of the load sharing 
function. The distributions have clearly two groups as 
found for other quantities also. 

In the limit where the long range interaction domi- 
nates, the clusters are randomly distributed on the lattice 
indicating that there is no correlated crack growth in the 
system as well as that the stress is not concentrated in re- 
gions. The cluster structure of the limiting case of 7 = 
can be mapped to percolation clusters on a square lat- 
tice generated with the probability < P(<J C ) < 1, where 
<7 C is the fracture strength of the fiber bundle. How- 
ever, the value of P(ct c ) depends on the Weibull index 
p and is normally different from the critical percolation 
probability p c — 0.592746 of the square lattice. Equal- 
ity P(a c ) — p c is obtained for p = 1.1132, hence, for 



physically relevant p values used in simulations the sys- 
tem is below p c at complete breakdown. This argument 
also justifies the exponential-like shape of the cluster size 
distributions of GLS in Fig. |[ This picture radically 
changes when the short range interaction prevails. In 
this case, the stress transfer is limited to a neighborhood 
of the failed elements and there appear regions where a 
few isolated cracks drive the rupture of the material by 
growth and coalescence. Thus, the probability of the ex- 
istence of a weak region somewhere in the system is high 
and a weak region in the bundle may be responsible for 
the failure of the material. The differences in the struc- 
ture of clusters also explain the lack of a critical strength 
when N goes to infinity in models with local rearrange- 
ment of stress. Since in the GLS model the clusters are 
randomly dispersed across the entire lattice, the system 
can "store" more damage or stress, whereas for LLS mod- 
els a small increment of the external field may provoke a 
run away event ending with the macroscopic breakdown 
of the material. 

Up to now, the change of the behavior of the system 
was observed for a certain value of 7 analysing various 
measured quantities. All these numerical results suggest 
that the crossover between the two regimes occurs in the 
vicinity of 7 = 2. Further support for the precise value of 
7 C can be obtained by studying the change in the cluster 
structure of broken fibers. The moments of n(s c ) defined 
as 

where is the fcth moment, describe much of the 
physics associated with the breakdown process. We will 
use these moments to quantitatively characterize the 
point where the crossover from mean field to short range 
behavior takes place. The zero moment n c = m is the 
total number of clusters in the system and is plotted in 
Fig. [ija as a function of the parameter 7. Figure fj]b rep- 
resents the variation of the total number of broken sites 
N c (the first moment N c — mi) when 7 increases. It 
turns out that up to a certain value of the effective range 
of interaction, N c remains constant and then it decreases 
fast until a second plateau seems to arise. Note that the 
constant value of N c for small 7 is in agreement with 
the value of the fraction of broken fibers just before the 
breakdown of the material in mean field models. This 
property clearly indicates a change in the evolution of 
the failure process and may serve as a criterion to calcu- 
late the crossover point. However, a more abrupt change 
is observed in the average cluster size (s c ) at varying 
7. According to the moments description, the average 
cluster size is equal to the second moment of the cluster 
distribution divided by the total number of broken sites, 
i.e. (s c ) = mijmx. It can be seen in Fig. 0d that (s c ) 
has a sharp maximum at 7 = 2.2 ± 0.1, and thus the av- 
erage cluster size drastically changes at this point, which 
again suggests the crossover point to be in the vicinity of 
7c = 2. 
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We now discuss the finite size scaling (FSS) of the 
avalanche distributions. For local load sharing one ex- 
pects that the cutoff in the avalanche distribution does 
not scale with L while for global load sharing the cutoff 
should scale with the size of the system. We have plot- 
ted in Fig. H the avalanche size distribution for several 
system sizes. As it can be observed, the FSS hypothesis 
is verified for the values of the exponent 7 correspond- 
ing to the global (Fig. |b) and the local (Fig. load 
sharing cases. Figure (Sc shows the moment analysis for 
five different system sizes in the range 2.0 < 7 < 2.5. 
It can be seen that the position of the maximum of the 
mi/mi curves is always at 7 = 2.2 ± 0.1, it does not 
scale with the system size. To determine the position 
of the crossover point more accurately we also analyzed 
the behaviour of a characterizing the strength of loga- 
rithmic size effect in Eq. (^), as a function of 7. From 
these studies it turned out that consistent interpretation 
of the numerical results can be given assuming that the 
crossover occurs in the vicinity of 7 C = 2.0 but stonger 
statement cannot be drawn due to the limited precision 
of calculations. 



IV. DISCUSSION AND CONCLUSIONS 



In the limiting case of global load sharing the break- 
ing of fibers is a completely random nucleation process, 
there is no correlated crack growth in the system, and the 
fiber failure which results in the catastrophic avalanche 
occurs at a random position in the system. As long as 
this microscopic damage mechanism holds when chang- 
ing the exponent 7, globally the system will behave in a 
global load sharing manner. On the other hand, when 
the load sharing is very localized, at the beginning of the 
failure process we get random nucleation of microcracks 
but later, correlated growth of clusters of broken fibers 
occurs. It then follows that along the perimeter of the 
clusters there is a high stress concentration and the final 
avalanche is driven by a fiber located at the perimeter 
of one of the clusters (the dominant one). At the fibers 
far away from the perimeter, the stress concentration is 
significantly lower, and the stress distribution is very in- 
homogeneous. In the case of localized load sharing this 
mechanism gives rise also to the logarithmic size effect as 
obtained also for the random fuse model Q . 

An interesting aspect to be explored in future work is 
whether there is a second transition point in the model 
when the 7-dependence of a c seems to disappear or it is 
just a crossover. In the case of localized load sharing the 
global failure is caused, as noted before, by the instabil- 
ity of a cluster of broken fibers, defining a critical cluster 
size. Furthermore, the exact value of 7 C might depend on 
the amount of disorder in the system which will also be 
subject of future studies. Preliminary studies of this issue 
indicate that the transition value 7 C gets slightly smaller 
as the system becomes more homogeneous (increasing p) . 



Before setting our conclusions, we would like to re- 
mark that a similar transition to the one obtained here 
has been observed in other models. It is well-known that 
long range interactions between spins can affect the crit- 
ical behavior of magnets. If we consider the addition of 
a term Ylrr' ~ r')s(r)s(r') to an Ising Hamiltonian, 
where V(r) ~ r ~( d + a ) corresponds to a long range inter- 
action decaying as a power law, a crossover to long range 
behavior is obtained provided that a < 2 — r/ sr , where 
rj sr is the value of the critical exponent 77 for short range 
force J39|-|d| . A similar behavior has also been reported 
i n p ercolation phenomena with long range correlations 

In summary, we have studied a fracture model of the 
fiber bundle type where the interaction among fibers is 
considered to decay as a power law of the distance from 
an intact element to the rupture point. Two very dif- 
ferent regimes are found as the exponent of the stress- 
transfer function varies and a crossover point is identified 
at 7 = 7 C . The strength of the material for 7 < 7 C does 
not depend on both the system size and 7 qualifying for 
mean-field behavior, whereas for the short range regime, 
the critical load vanishes in the thermodynamic limit. 
The behavior of the model at both sides of the crossover 
point was numerically studied by recording the avalanche 
and the cluster size distributions. The numerical results 
suggest that the crossover point falls in the vicinity of 
7 C = 2.0. Finally, we have outlined some general ideas 
which will be the subject of a forthcoming publication. 
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FIG. 1. Illustration of the model construction. X in- 
dicates a fiber, which is going to break, and O is an intact 
fiber in the square lattice. 
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FIG. 2. Ultimate strength of the material for different sys- 
tem sizes as a function of the effective range of interaction 7. 
A crossover from mean field to short range behavior is clearly 
observed. 
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FIG. 3. Variation of the material strength with N for sev- 

eral values of 7 . Note that when 7 increases the critical load plG g Snapshots of the clusters just before the comp i et e 

vanishes in the thermodynamic limit, whereas, for small 7 it breakdown of the mate rial. The change in the structure of 

has a nonzero value independent on the system size. tfafi dusterg can bg geen The vames of ^ &re . a ; 7 = 0] b) 

7 C = 2.2, c) 7 = 3, and d)-y = 9. 
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FIG. 4. Avalanche size distributions for different values of 
the exponent of the stress-transfer function 7. The upper 
group of curves can be very well fitted with a straight line 
with a slope r = — | (L = 257). 



FIG. 6. Cluster size distributions for different values of 
the stress-transfer function exponent 7. Clearly, two differ- 
ent groups of curves can be distinguished as found for other 
quantities also (L = 257). 
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FIG. 7. Moments of the cluster size distribution as a func- 
tion of 7 (see text for details on the definition of ink)- A sharp 
maximum is observed at 7 = "/ c ~ 2.2 for the average cluster 
size (s c ) = =2. 
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FIG. 8. Finite size scaling analysis, a) Scaling of the cutoff 
with the system size for the local load sharing case, b) Scaling 
of the cutoff with the system size for the global load sharing 
case, and c) Average cluster size, (s c ) = 2^-, for different sys- 
tem sizes. Note that in c) the position of y c does not change. 
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